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Abstract 

The  limitation  of  the  speed  and  memory  of  calculating 
machines  places  an  upper  bound  on  the  number  of  meshpoints 
that  may  be  used  in  a  finite  difference  calculation.   This 
means  that  in  problems  involving  many  independent  variables 
(and  for  present-day  machines,  three  is  many)  the  mesh 
employed  is  necessarily  coarse.   Therefore  in  order  to  get 
reasonably  accurate  final  results  one  must  employ  highly 
accurate  difference  approximations .-   The  purpose  of  this 
paper  is  to  set  up  and  analyse  such  difference  schemes  for 
solving  the  initial  value  problem  for  first  order  symmetric 
hyperbolic  systems  of  partial  differential  equations  in  two 
space  variable. 
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DIFFERENCE  SCHEMES  WITH  HIGH  ORDER  OP  ACCURACY  FOR 
SOLVING  HYPERBOLIC  EQUATIONS. 

Introduction. 

It  Is  well  known  that  a  difference  scheme,  In  order  to 
furnish  accurate  answers  over  a  reasonably  long  range  of 
time  has  to  be  stable.   The  bulk  of  this  paper  Is  devoted 
to  determining  the  conditions  under  which  the  proposed  dif- 
ference schemes  are  stable,, 

In  Section  1  we  give  a  brief  review  of  the  general 
theory  of  accuracy  and  stability  of  difference  schemes.   In 
Section  2  we  show  that  the  set  of  matrices  whose  field  of 
value  belongs  to  the  unit  circle  forms  a  stable  family. 
This  result  Is  of  Interest  In  Its  own  right.   In  Section  3 
we  set  up  some  difference  schemes  of  second  order  accuracy 
and,  with  the  aid  of  the  criterion  described  In  Section  2, 
analyse  the  range  of  parameters  for  which  these  schemes  are 
stable.   In  Section  4  we  give  a  geometric  Interpretation  of 
the  stable  range  of  the  parameters  and  In  this  connection 
we  devise  a  difference  scheme  with  a  maximum  stable  range; 
this  scheme  however  Is  accurate  only  to  first  order.   Section 
5  contains  some  remarks  and  open  questions  concerning  the 
effect  on  stability  of  the  non-constancy  of  the  coefficients 
of  a  difference  scheme.   In  Section  6  we  show  how  to  set  up 
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difference  schemes  with  higher  order  accuracy  for  nonlinear 
hyperbolic  systems  of  conservation  laws,  such  as  the  equations 
of  compressible  flow  and  magnetohydrodynamlcs. 

The  difference  schemes  discussed  In  this  paper  are  at 
present  being  tried  out  on  various  problems  by  colleagues  at 
the  AEC  Computing  Center  at  N.Y.U.  and  at  the  Los  Alamos 
Scientific  Laboratory.   We  hope   that  the  reports  of  these 
practical  evaluations  of  the  methods  proposed  here  will  soon 
be  available. 


•X- 

We  also  hope  that  they  will  indicate  that  our  methods  are 
useful  in  practice. 
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1.   Review  of  the  Notion  of  Stability  and  Accuracy. 

The  class  of  equations  under  consideration  are  of  the  form 


( 1 o 1 )        u ,  =  Au  +  Bu  , 
^    ^         t     X     y' 


u  a  vector  function  of  x,  y  and  t   and  A   and  B   symmetric 
matrices  which  may  depend  on  x   and   y|  for  sake  of  convenience 
we  shall  not  consider  explicit  dependence  on  t .   On  occasion  we 
shall  abbreviate  the  right  side  of  ( 1.1 )  by   G  and  write  the 
equation  In  the  form 


(l.l)        u,  =  Gu^ 


and  Indicate  explicitly  only  the  dependence  of   u   on   t. 

We  are  Interested  In  the  Initial  value  problem^  I.e.  of  finding 

a  solution  of  (l,l),  given  the  value  of   u(o). 

We  shall  consider  difference  approximations  to  (l.l)  of  the 
form 


1.2)        v(t  +  h)  =  S^v(t); 


here   v   denotes  an  approximation  to   u,   h   Is  the  time 
Increment  and   S   Is  a  difference  operator: 


(lo)         S,  =  V  c  .  t\ 
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where   j   Is  a  multiindex   (j"  ,j\  )   and   T'^   abbreviates 

■'^    y 

T    T  /   where   T   and   T   denote  translations  by  the  amounts 
-^    y         X      y 

\xh     and   vh   In  the   x  and  y  directions  respectively,  \i 

and   V   being  constants  Independent  of  h  .   The  coefficient 

matrices   c.   are  functions  of  x,y  and  of  h;  they  are  polynomials 

In  h  . 

Definition  :  (I.2),  (lo)   approximates  ( 1 , 1 )  with  m 

order  accuracy  If  for  all  smooth  solutions   u(t)   of  (l.l) 


(1.^)        f|u(t  +  h)  -  S^u(t)![  <  o(h"^+i). 

I.e.    If  exact  and  approximate  solutions  differ  after  one  time 
step  only  by   0 ( h'^^  "  ) . 

Definition :   (I.2)   Is  stable  If  Its  solutions  are  uniformly 
bounded  In  a  unit  time  range,  I.e.  If  there  exists  a  constant   K 
such  that 


s[;f  f  <  K 


for  all  n,h   satisfying  nh  <  1  . 

In  this  paper  the  norm  appearing  in  these  definitions  will 
be  taken  as  the  L   norm. 


The  following  is  well  known,  see  [l^^]: 
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Theorem  1.   Let   u   and   v   denote  solutions  of  the  exact 
differential  equation  (l.l)  and  the  difference  equation  (1.2) 
respectively,  having  the  same  smooth  Initial  values.   Then: 

(1.5)      I [u(t)  -  v(t)| I  <  0(h™),   t  <  1  . 

for  all  smooth  Initial  values  If  the  difference  scheme  Is  stable 
and  accurate  of  order  m. 

Since  accuracy  and  some  sort  of  stability  are  necessary  as 
well  as  sufficient  In  order  that  the  overall  error  be  of  the 
order  (I.5),  we  endeavor  to  construct  difference  approximations 
accurate  of  order  m  and  also  stable.   In  the  present  paper  we 
take   m   to  be  2. 

It  Is  easy  to  express  accuracy  of  order  m  as  a  certain 
number  of  linear  relations  among  the  coefficients   c.   of   S 
(see  [4]).   By  Taylor's  theorem 


,  m  .  T 

u(h)  =  u(o)  +  hD,u(o)  +    ...    +  ^   D%(o)  +  0(h'^+^) 


where  D  denotes  differentiation  with  respect  to  t.  The  dif- 
ferential equation  (l.l)'  asserts  that  for  solutions  of  (l.l)', 
D   =  G;   so  for  solutions   u: 

m 

(1»6)     u(h)  =  YZ     ^  gMo)   +  0(h"^+^). 

n=o 


On  the  other  hand  expanding  T'^v   in  a  Taylor  series  around 
(x,y)  gives 


T^'v  =  M.v  +  0(h"^'*"^) 


where   M .   is  a  partial  differential  operator  of  order  m  whose 
coefficients  are  polynomials  in  h   of   order  m.   Substituting 
into  (1.3)  we  get 

(1.7)        v(h)  =  Mu  +  0(h"'+^), 

where   M   is  a  differential  operator  as  above  whose  coefficients 
depend  linearly  on  the   c  .,   Equating  the  coefficients  of  M  with 
those  of  the  operator  on  the  right  of  (I.6)  gives  the  consistency 
conditions   which  are  necessary  as  well  as  sufficient  for  m 
order  accuracy. 

There  is  a  more  transparent  form  of  writing  these  consistency 

conditions;  we  shall  derive  this  form  first  in  the  case  when  the 
coefficients  of  the  differential  and  difference  equations  are 
constant,  i.e.  independent  of  x   and   y  .   We  start  with  the 
observation  that  it  suffices  to  verify  the  error  estimate  (1.5) 
for  a  dense  set  of  solutions.   We  choose  these  solutions  as  the 
exponential  ones;  that  is,  prescribe  u(o)   as 


u( 


,(0)  ^  e^^^^+y^V 

where  ^,x\      are  arbitrary  real  numbers  and  ^     is  an  arbitrary 
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vector.   The  corresponding  solution  of  (l.l)  Is,  it  is  easy  to 
check 


(1.8)        u(t)  =  e"^^^+^'^)u(o 


while  the  corresponding  solution  of  (1.2)  is 


(1.9)        v(h)  -  C((i.h4,vhr|)u(o) 


where 


(1.10)  C(C)  =2_  c.e^J^. 

Here   C   denotes  a  vector. 

Comparing  (I.8)  and  (1.9)  we  see  that   accuracy  of  order 
m  means  that 

(1.11)  e^(^^+^l^)  .  C(^^,v,)  -  0(lcr'+^) 

for   C   near  zero. 

For  equations  with  variable  coefficients  it  is  not  hard  to 
show  that  for  accuracy  of  order  m  we  must  have 

(1»11)'       e^(^^+^^^  =   C^(ne,vrJ  +  Oil^ir"-^ 


for  every  x,      where   C   =  C  (x,y,0   is  defined  as  the  value 
of   C   at   h  -  0  . 
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The  function   C   defined  by  (l.lO)  is  called  the  amplifi- 
cation matrix  of  the  difference  operator  (lo2). 

We  turn  now  to  the  question  of  stability  of  difference 
schemes >   We  shall  deal  here  with  the  case  of  constant  coef- 
ficients; schemes  with  variable  coefficients  will  be  discussed 
in  Section  5. 

Denote  the  Fourier  transformation  in  the  space  variables 
by  T;    then 

TS^v  =  C(hC)Tv, 

where   C   denotes  the  dual  variable.    Repeated  applications  of 
the  above  identity  gives 

TS^v  =   c"(hC)Tv. 

Since  Fourier  transformation  is  an  isorrietry  In  the  L   norm, 

n 
the  uniform,  boundedness  of   S,  v   is  equivalent  to  the  uniform 

boundedness  of  their  Fourier  transforms.   The  latter  is  clearly 

equivalent  to  the  uniform  boundedness  of  the  matrices   c".   Thus 

we  have  shown : 

Theorem  2.   A  difference  scheme  with  constant  coefficients 
is  stable  if  and  only  if  all  powers  of  the  associated  amplification 
matrix  are  bounded,  uniformly  for  all  real  values  of  C   and  all 
powers  of  the  matrix. 

To  be  able  to  use  the  above  stability  criterion  we  need  to 
know  conditions  under  which  a  family  of  matrices  has  the 
property  that  all  powers  of  its  elements  are  uniformly  bounded. 
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As  observed  by  von  Neumann,  a  necessary  condition  for  this  Is 
that  the  eigenvalues  of  each  matrix  of  the  family  be  not  greater 
than  one  In  absolute  value.   This  condition  by  Itself  is  not 
sufficient;  there  are  various  additional  conditions  given  in 
the  literature,  see  [  1^  ]  and  [  5  ]  which  together  with  von 
Neumann's  condition  guarantee  the  uniform  boundedness  of  the 
set   <:  c"  f  ,   Necessary  and  sufficient  conditions  were  given  by 

L    J' 

Kreiss  [  6  ],     [  7  ]  and  by  Buchanan  [  2  ],   In  the  next  section 
we  shall  give  a  new  sufficient  condition  and  use  it  In  Section 
3  to  discuss  the  stability  of  the  difference  schemes  which  are 
the  subject  of  this  paper. 
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2 .    A^  Stability  Theorem . 

Theorem  3.   Suppose  that  the  field  of  values  of  a  matrix 
C   lies  in  the  unit  disk,  i.e.  that 

(2,1)         I  (Cu,u)|  <  1 

for  all  unit  vectors   u.   Then  there  exists  a  constant   K 
depending  only  on  the  order  of   C   such  that 


C"|  <  K,         n  =  1,2 


Here   fc  f   denotes  the  operator  norm  of   C   with  respect 
to  the  Euclidean  norm  for  vectors. 

Rem.ark :   Since  all  eigenvalues  of   C   belong  to  its  field  of 
values,  (2.1)  Implies  that   C   satisfies  the  von  Neumann  condition, 

Proof  of  theorem  3=   It  is  well  known  that  every  matrix  is 
unitarlly  equivalent  with  an  upper  triangular  matrix.   Since  the 
field  of  values  of  two  matrices  so  related  is  the  same,  and  since 
the  norms  of  their  powers  are  the  same,  we  may  assume  that   C 
is  already  upper  triangular. 

We  shall  treat  first  the   2  x^  2  case: 


a   b  \ 

c  =        : 

0   c  / 
\ 

Then 

(Cu,u)  =  a[t[   +  b  t  z  +  c[z 
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where   t,z   are  the  components  of   u  .   Choose   |t|   =  |zi   =2 
and  adjust  the  argument  of   t   and   z   so  that  the  middle  term 
has  the  same  argument  as  the  sum  of  the  first  and  the  third. 
With  such  a  choice 


n         \ I    I a+C|  ,   b 
Gu,  u  j  I  =  I  -—5—1  + 


2  '     2 
So  from  (2«l)  we  have 

(2.2)  I b|  <  2  -  I a+c|  , 

By  the  triangle  inequality   | a+c [  >  2|c[  -  [a-ch  thic  and  (2,2) 
imply  that 

(2.3)  |t)|  <  2(  [1  -  I  c|  )  +  [a-c|  . 

Consider  now  powers  of   C;  it  is  easy  to  show  recursively 
that 

a     P  (a,c  )lD 
{2A)  C"  = 

0 


^    I         \  n-1    n-2  n-1   a  -c 

where   P  =  P  (a,c)  =a    +ac+...+c    =  ^  _  ^ 


As  remarked  before,  it  follows  from  (2ol)  that  the  eigenvalues 
a   and   c   of   C   do  not  exceed  one  in  absolute  value.   So 


-  V\ 


using   |a|  <  1   we  have  from  the  first  form  for   P  the  inequality 


(2.5)      |P|<l+|c|  +  ...  +  |cr^^<    ^ 


1-  c 


Prom  the  second  form,  using   |a|  ,   |c|   <  1  ,  we  get 


(2.6)        |P|  < 


2 


a-c 


Multiplying  (2.5)  by   2(l  -  |c|),   (2.6)  by   |a-c|   and  adding 
the  two  we  get 

(2.7)  [2(1  -  |cl )  +  |a-c| ]  lP|  <  4  . 
Combining  (2.3)  and  (2.7)  we  get 

(2.8)  |bP|  <  4 

n 
which  shows  that  the  corner  element  of   C   is  at  most  four. 

The  above  derivation  of  (2.8)  from  (2.3)  is  taken  from 
Buchanan ' s  paper  [  2  ] . 

We  shall  now  prove  theorem  3  for   pxp  matrices  C 
inductively  on  p;  this  device  is  borrowed  from  de  Bruijn  who  has 
used  it  in  [  1  ] . 

Let   C  be  an  upper  triangular   pXp  matrix.   We  write  it 
in  block  notation  as 


15 


A   p 
o   c 


where  A   is   a(p-l)  X  (p-l)   upper  triangular  matrix,   p   is 
a  column  vector  with   p-l   components,  and   c   is  a  scalar. 

Let   u   be  a  column  vector  with  p   components  and  length 
1  .   We  can  write  it  as 

■  t  V  ■' 
u  = 

where   v   is  a  unit  vector  with   p-l   components  and   t,z  are 
complex  numbers  satisfying 


tl^  +  Izl^  =  1  . 


With  this  notation,  the  field  of  values  of  our  matrix   C   may 
be  written  in  the  form 


(Cu,u)  =  (Av,v)|t|^  +  (p,v)tz  +  c|z|^ 


=   a  I  tl  ^  +  b  tz  +  cl  zl  ^ 


where   a   and   b   abbreviate 


(2.9)         a  =  (Av,v),       b  =  (p,v) 
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since   C   satisfies  (2.l)  It  follows  that  the  absolute  value  of 
the  expression   a|t|   +  "bt  z  +  c[z|    does  not  exceed  1.   Since 
this  expression  can  be  thought  of  also  as  the  field  of  values 
of  the   2\2  matrix 


we  conclude  that   a,b,c   satisfy  the  Inequality  (2o3)  where  we 
now  substitute  the  expressions  (2.9)  for   a  and   b  <,   This 
yields 


(2,10) 


(P,v)l  <  2(1  -  |cl)  +  |((A-cl)v,v) 


here  we  wrote   a-c  =•  (Av,v)  -  c  =  (Av,v)  -  c(v,v)  =  ((A-cI)v,v), 
since   V   is  a  unit  vector.   We  also  conclude  that   |a|  <  1, 
i.e.  that  the   (p-l)  X  (p-l)   matrix  A   satisfies  the  hypothesis 
(2ol).   Moreover,   |c|   does  not  exceed  1. 

Our  aim  is  to  derive  a  uniform  bound  for  all  powers  of   C. 
It  is  easy  to  verify  recursively  that  these  may  be  written 


■.n 


P  (A,c)p 


n 


where   P  (A.c)   is  the  matrix 
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(2.11)  P  =  a"-^  +  ca"-2  +  ...  +  c"'4  -  (a"-c"i)(A-cI)"\ 

By  the  induction  hypothesis  there  exists  a  constant   K   such 
that 

(2.12)  I  A^^l  <  K  ,      n  -  1,2,  .  .  . 

t  n  I 
Thus  in  order  to  find  an  upper  bound  for   | C  |   it  suffices  to 

find  one  foi'  the  norm  of  the  vector   Pp . 

Prom  the  first  expression  for   P   in  (2.11 ),  using  (2.12), 

we  find  that 


(2.13)    fP|  <  [A*^-^!  +  ...  +  Ic""-"^!  <  K(l  +    ...    +    IcT-)   ^ 


1-  c 


Next  we  note  that  according  to  a  well-known  principle  the 
norm  of  the  vector  Pp  can  be  characterized  in  terms  of  inner 
products : 

(2ol4)   I  Pp  I  -  Sup    |(Pp,w)!  =   Sup   I  (p,P'''w)[  , 
[  w  I  =  1  I  w  I  =  1 

Thus  we  can  bound   | Pp [   from  above  by  finding  an  upper  bound 
for   (Pp,w)  =  ((3,P-^-w). 

We  now  choose  the  unit  vector 

V  =^  P*w/|  P-'-w|  ; 


then   (2.10)  becomes 

(2.15)   I  (p,P*w)|/|P*wI  <  2(l-|c|)  +  I  ((A-cl)P*w,P*w)|/|P*- 

The  second  term  on  the  right  can  be  rewritten  as 

(P(A-cl)P*w,w)/| P*w| ^; 

and  using  the  second  expression  In  (2.II)  for   P  we  can 
rewrite  this  as 


w 


(2.16)         ((a"-  c"i)P*w,w)/Ip* 


w 


By  (2.12)  and  since   Ic|  <  1   the  norm  of  the  operator 


A   -  c  I   Is  at  most   K  +  Ij   so 


I  I  (a" -c"l)P*w[  I  <  (K+l)|  |P*w|  I  . 

Therefore  if  we  estimate  the  numerator  in  (2.I6)  by  the 
Schwarz  inequality^  then  using  the  above  estimate  and  the 
fact  that   w   is  a  unit  vector  we  find  that 


K+1 


P*w 


is  an  upper  bound  for  (2.I6).   Substituting  this  upper  bound 
for  the  second  term  on  the  right  in  (2.I5)  we  get,  after 
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multiplying   by   |  P'^'w|  : 

I  ( 3 , P*w ) [  <  2 ( 1- I c I  )  I P*w I  +  K  +  1 . 

Since   w   is  a  unit  vector   | P*wl  <  Ip*|,  the  norm  of   P* 
equals  that  of   P,   for  which  we  already  have  the  estimate 
(2.13)j  thus  we  obtain 

I (p,P*w)[  <  3K  +  1    for  all  unit  vectors   w. 

In  view  of  (2.14)   the  above  inequality  shows  that 


P^PI  <  3K  +  1; 


this  is  the  required  uniform  bound,  and  the  induction  is 
now  completed. 

Observe  that  the  value  of  the  uniform  bound   K 
derived  here  depends  on  the  order   p   of  the  matrix   C, 
and  increases  exponentially  with   p  .   It  would  be  of  some 
interest  to  study  the  dependence  on   p   of  the  best 
constant   K   in  theorem  3- 

Corollary  to  theorem  3:   Let   C   be  a  matrix  satisfying 
the  hypothesis  of  theorem  1,  and  let   A   be  any  matrix 
whose  norm  does  not  exceed  M,   Then  for  every  positive   h 


(2.17)         I  (C  +  hA)"|  <  K(l  +  hM)", 
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To  see  this  we  merely  observe  that  the  matrix 

(C  +  hA)/(l  +  hM)   also  satisfies  the  hypothesis  of  theorem 

1. 

The  Inequality  contained  In  this  corollary  can  be 
used  to  study  the  stability  of  difference  schemes  for 
Inhomogeneous  differential  equations  with  constant   ^  '..  _ 

coefficients,  and  also  the  stability  of  difference  schemes 
with  variable  coefficients. 
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3  .    Derivation  and  Analysis  of  Difference  Schemes  _  o  _f 

Se^ond^Jlrder  _A  ccur  a  c  y . 

We  shall  construct  now  with  the  aid  of  condition  (l.ll) 
some  difference  schemes  which  are  accurate  to  second 
order.   We  take   A   and   B  to  be  constants  and  for 
simplicity   n  -  V  =  1.   Expanding  In  Taylor  series  near 
^  =  ri  =  0  we  have 


(5.1)    e^^^^+^^^  =  I  +  l(^A+i-iB)-  i  (lA+TiB) 


2 


where  the  symbol  ~     denotes  congruence  modulo  third  order 
terms.   Furthermore 

t  —  sin  I,   T]  =  sin  1-1  ,  iri   —sin  i    sin  ^ 
(5.2) 

^2/2  -  1  -  cos  i,      T]  /2=  1  -  cos  T]. 

Substitute  the  congruences  (3.2)  Into  the  right  side  of 
(3.1)  and  denote  the  resulting  function  by   C(|,ti): 

(5^5)    c  =  I  +  l(A  sln^  +  B  sin  ii )  -  k^{l    -    cos  i) 
-  i  (AB  +  BA)  sin  |  sin  t]  -  B^(l  -  cos  r, ) . 

By  the  above  construction  and  (l.ll)   C   Is  the  amplifi- 
cation matrix  of  a  difference  scheme  which  Is  accurate  to 
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second  order. 

Clearly   C   as  given  by  (3o3)  is  a  trigonometric 
polynomial  of  degree  one  in  each  variable.   This  means 
that  If  we  write   C   In  the  form  (l.lO),  the  range  of 
the  multllndex   j   Is   -  1  <  j  ,j   <  1  .   Thus  the 

^  y 

operator   S   in  (lo)  can  be  expressed  in  terms  of 
translations  by  the  amount   ±  h   in  the  x  direction 
and   t  h   in  the   y  direction »   This  means  that  the 
difference  equation  (l.2)  expresses  the  value  of   v   at 
the  point   (x,y)   and  time   t   in  terms  of  values  of   v 
at  time   t   and  the  points  {x,j)      and  their  eight  lattice 
neighbors,  as  shown  in  figure  1: 


Figure  1 

For  this  reason  the  difference  scheme  associated  with 
(5.5)  is  called  a  nine-point  scheme. 

(5.3)  is  not  the  only  nine-point  scheme  which  is 
accurate  to  second  order. 

The  scheme  associated  with  the  amplification  matrix 


(30) 


C   =  C  - 


A^+B^ 


(1  -  cos  ^)(l  -  cos  T])  , 
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C   as  given  by  (3o)  has  the  same  accuracy  as   C,  since 
the  added  term  In   C    Is  of  fourth  order. 

Theorem  4 :   The  difference  scheme  associated  with 
(3-3)  is  stable  if 


(3.4)  a2  <  1  I,    b2  <  1  I 


and  the  scheme  associated  with  (3=3)    is  stable  if 


(3.4)'  A^  +  B^  <  i  I 


These  results  are  the  best  possible  ones  in  the  sense 
explained  at  the  end  of  the  proof. 

Proof:   We  shall  show  that  under  the  conditions  stated 
above  both   C   and   C    satisfy  the  hypothesis  of  Theorem  3 
We  start  with   C;  separating  it  into  real  and  imaginary 
part  we  write 

(3.6)  C  =  R  +  iJ, 

Here 

(3.7)  J  =  A  sin  4  +  B  sin  r\ 

and 

(3o8)         R  =  I  -  K, 
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where   K   abbreviates 


.2 


(3.9)  K  =  A^(l-  cos  i)   +   B^(l-  cos  ri)  +  ^(AB+BA)sln  ^  sin  t^  , 

It  turns  out,  and  this  Is   of  Importance,  that   K  can  be 
written  as  a  sum  of  three  squares.   Using  the  abbreviations 

(5.10)  1  -  cos  I  =  X     ,      1  -  cos  Tj  =  Y 

we  have  the  Identity 

(5.11)  K  =  I  A^X^  +  i  B^Y^  +  i  J^  . 

The  verification  of  this  identity  is  left  to  the  reader. 

Clearly  since   A   and  B  are  symmetric  so  are   R 
and   J. 

Our  aim  is  to  estimate  the  quantity  (Cu,u)   for  unit 
vectors   u  .   We  can  write 

(Cu,u)  =  (Ru,u)  +  i('ju,u)  =  r  +  Ij.- 

Since   R  and   J  are  real  and  symmetric,  r  and   j   are 
real,  and  so 


(3.12)  f(Cu,u)|2  =  r2  +  j2 


We  estimate   j   by  the  Schwarz  inequality: 
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(3.13)         j^  =  (Ju,u)^  <  I  Jul^ 


By  (3-8)  we  can  write 
(3.1^)  r   -  (Ru,u)  =  1  -  (Ku,u) 

By  (3.11)  we  can  write 

(3.15)  (Ku,u)  =    i  |Au|^X^  +  i  |Bu|^Y^  +  ^  \Ju\^ 


i  a  X^  +  I  b  Y^  +  i  I  Jul 


where  we  have  used  the  abbreviations 


(3-16)         |Au|  ^  =  a,        |Bu|  ^  =  b. 


Squaring  (3.1'^)  and  using  (3. 15)  we  get 


r^  =  1  -  a  X^  -  b  Y^  -  |Ju|^  +  (Ku,u)^. 


Adding  (3. 13)  to  this  we  get 


(3-17)     r^  +  j^  <  1  -  a  X^  -  b  Y^  +  (Ku,u)^, 


Next  we  turn  to  estimating   (Ku,u);  using  (3-9)  we  have 
(Ku,u  )  =  aX  +  b  Y  +  2Re(Au,Bu)sin  ^  sin  r\. 
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Applying  the  Schwarz  Inequality  to  the  last  term  on  the 
right  we  get 

(3  =  18)  I  (  Ku ,  u )  [  <  a  X  +  b  Y  +  |  Au  |  |  Bu  1  |  sin  1 1  [  sin  r]  |  . 

Estimating  the  last  term  in  the  above  by 

|Au|^sin^4  +  iBul^sin'^Ti 

and  using  the  elementary  inequalities 

•  2* 
sir.  £  ^  -,       f      ^ 
— 7:r^   <  1  -  cos  t,  etc. 
d.        — 

we  get 

1 (Ku,u)|  <   2aX  +  2bY. 
So 
(3.19)    (Ku,u)2  <  Sa^X^  +  Sb^Y^; 

substituting  this  into  (3.I7)  gives 

(5-20)    r^  +  j^  <  1  -  aX^d  -  8a)  -  aY^(l  -  8b). 

The  expression  on  the  right  will  not  exceed  one  if  a   and 
b  are  both  not  greater  than  -q  .   According  to  the 
definition  of   a  and   b   this  is  the  case  If  condition 
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(3o4)   of  theorem  4  is  satisfied.   Thus  combining  (5-19) 
and  (3.12)  we  see  that  under  condition  (3-^)  the  field 
of  values  of   C   lies  in  the  unit  disk.   According  to 
theorems  2  and  5  this  proves  the  stahlllty  of  the 
associated  scheme.. 

Observe  that  in  the  estimates  above  the  Schwarz 
inequality  was  used  so  gently  that  the  sign  of  equality 
can  hold  throughout.   In  fact  if   A  =  B  one  can  easily 
show,  by  setting  i   =   n    ,    that  condition  (3.4)  is 
necessary  as  well;  in  this  sense  our  result  is  best  possible 

We  turn  now  to  the  second  part  of  theorem  4. 
Following  the  same  line  of  argument  we  get  in  the  place  of 
(3.9)  and  (3.11); 


k'  =  K  +  ^  XY. 


In  place  of  (3-1?)  we  get 
(3.17)'   ^2  ^  j2  ^  ^  _  g_x2_  bY^  -  (a+b)XY  +  (Ku,u)' 


=  1  -  (aX  +  bY)  (X+Y)  +  (Ku,u)^ 


and  in  place  of  (3.I8)  we  get,  after  replacing   |AuI|Eu| 

■K    a+b  . 
by  -g-  • 

(3.18)'     1  (Ku,u)|    <    aX  +   bY   +  ^      ■,  |sin    ll  |sin   ril    +   XYC 
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We  estimate  the  curly  brackets  by  the  Schwarz  Inequality; 


<  ^jsln^|  +  (1-  cos  4)'  ^ 


f2(l-  cos4)  ~S 


=  2  ^rxY 


Substituting  this  into  (3«l8)   gives 
(3.18)"    I  (Ku,u)|  <  (a  VX  +  b/Y)(^X  +  ^). 
Squaring  and  using  the  Schwarz  inequality  gives 


(3.19)'    (Ku,u)^  <  2(aX  +  bY)(a  +  b)(X  +  Y), 


Substituting  this  into  (3«17)   gives 

(3.20)'    r^  +  j^  <  1  -  (aX  +  bY)(X  +  Y)[l  -  2(a+b)]. 

Clearly  the  expression  on  the  right  will  not  exceed  one  if 
2(a  +  b)  does  not  exceed  one;  but  this  is  precisely  what  is 
guaranteed  by  condition  (3-^)  .   This  completes  the  proof 
of  the      d  half  of  theorem  4. 

Again,  by  setting  A  =  B  and   ^  =  ri   one  can  easily 
show  that  condition  (3-4)  also  is  necessary. 

So  far  •  fj.  =  V  =  1;   it  is  easy  to  show 

that  in  g'   .  .1  (3.''0  and  (3. 4)   have  to  be  formulated  as 
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follows 


a2    IB-    I     ,    A^  B^  ^  I 

Observe  that  condition  (5.^)  'is  less  restrictive  than 
(5.4),  i.e.  that  the  scheme  associated  with   C    is  stable 
in  a  wider  range  than  that  associated  with   C.   This 
greater  stability  is  the  effect  of  the  extra  term  in   C 
which  Introduces  in  the  associated  difference  operator  a 
corresponding  additional  term.   This  additional  term  is 
clearly  a  difference  analogue  of  the  operator 

_  h^  a!+b!  ^2  ^2  ^ 

2     X  y 

Such  a  hiR.her  order  negative  definite  term  is   ;.:. 
called  an  artificial  viscosity.   The  effect  of  such  terms 
has  been  investigated  by  several  authors,  in  particular  by 
Kreiss  [  7  ]   in  the  linear  case  and  by  von  Neumann  and 
Rlchtmyer  [  13 ]   and  Lax  and  Wendroff  [  11 ]   in  the  non- 
linear case.   The  results  of  this  section  furnish  another 
Illustration  of  the  stabilizing  effect  of  artificial 
viscosity. 
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h.        The  Geometric  ffeanlng  of  Stability. 

The  difference  schemes  discussed  in  the  last  section 
are  only  conditionally  stable,  i.e.  they  are  stable  only 
if  the  coefficients  of  the  differential  equation  which 
they  approximate  satisfy  the  inequalities  {^.^)    and  (5.4)  . 
An  intuitive  reason  why  such  inequalities  are  necessary  for 
stability  has  been  given  a  long  time  ago  by  Courant, 
Friedrichs  and  Lewy  in  their  classical  paper  :      ; 

Let   p  be  any  point  and   t   any  time;  denote  by 
D(p,t  )  the  set  of  those  points  on  the  Initial  plane   t  =  0 
where  the  values  of  the  initial  data  influence  the  value  of 

the  solution  of  the  differential  equation  at   P,t   . 

o 

Denote  by  D,  (p,t  )   the  analogous  set  with  respect  to 

the  difference  equation.   Then,  as  Courant,  Friedrichs  and 

Lewy  have  pointed  out,  if  a  difference  scheme  is  convergent 

for  all  smooth  initial  data  then  for  any   p  and   t    the 

set   D(p,t  )   must  be  contained  in  the  set  of  limit 

points  of  D  (p,t  )   as   h   tends  to  zero.   Since  convergence 

and  stability  are  equivalent  (theorem  l)  this  gives  a 

necessary  condition  for  stability. 

To  use  this  condition  we  have  to  determine  the  domains 

of   dependence    D  and   D,  .   If  we  deal  with  the  case 

h 

of  constant  coefficients  then  by  reason  of  homogeneity  we 
may  take  p  to  be  the  origin  and   t    to  be  one.   Taking 
(i,  =  V  =  1,   for  a  nine-point  scheme  the  set   D   consists 
then  of  all  lattice  points  with  mesh  width  h  inside  the 
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unit  square  S. 


A 


y 


/ 


X 


Figure  2 

The  set  of  limit  points   of  D   is  therefore  the  unit  square   S, 
The  determination  of  the  domain  of  dependence   D   Itself  Is 
a  slightly  delicate  problem;  but  the  support  function  of   D 
is  easily  determined.   We  recall  that  the  support  function 
hpj  ( ^ ,  1] )   of  any  closed  bounded  set  in  the  plane  is  defined  as 
follows : 


hj^(^,n)  =  Max        (x^  +  yrj) 
(x,y)  in  D 


A  well  known  result  in  the  theory  of  hyperbolic  equations 
(see  e.g.  [  9  ])  asserts: 


(4.1)     h^(^,ri)  =  A    (^A  +  iiB) 


Where  '^^qjX^)      denotes  the  largest  eigenvalue  of  X. 

(This  result  is  related  to  the  fact  that   A  ^^(^A+i-iB)   is 
the  majclmum  speed  of  propagation  in  the  direction   (I,t])). 
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It  follows  from  the  definition  of  support  function 
that  if  one  set  is  contained  in  a  second  set,  the 
support  function  of  the  first  does  not  exceed  that  of 
the  second.   It  is  further  known  from  the  theory  of 
convex  sets  that  if  the  second  set  is  convex,  then  the 
converse  of  the  above  statement  holds.   Thus  we  can 
express  the   C-P-L   condition  in  this  form: 

A  nine  point  scheme  for  equation  (l.l)  can  be  stable 
only  if 


(^.2)  %  <'^S 


for  all  i,r\,      where   S   denotes  the  unit  square. 

It  is  easy  to  show  that  if  (4.2)   holds  for  all 
(^,'ri)   which  are  perpendicular  to  a  side  of   S   then 
it  holds  for  all  ^,r\.      Thus  taking   (^,rj)   to  be 
(^  1,0)   and   (0,^l)  we  see,  using  (4,2),  that  the 
necessary  condition  of   C-P-L   can  be  stated  as  follows 


1  <  >  .   (A)  <  A    (A)  <  1 
—  mm    ^  —     max^  '  — 


1  <  A  .  (3)  <  A     (E)  <  1, 
—  mm^  ^   —  max  ^  ^  — 


These  Inequalities  about  eigenvalues  can  be  expressed 
also  as  matrix  '    'ilitles: 

(4.3)  A^  ^  I  ,    B^  <  I 
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It  is  something  of  a  curiosity  that  both  {^.h)    and  (3.'^)   are 

more  stringent  than  (4.3). 

It  is  easy  to  give  a  geometric  interpretation  of  (3- 4)  and 

(3- 4)  .   The  former  asserts  that  the  domain  of  dependence   D 

t 
lies  inside  the  small  square  shown  on  Figure  3j  while  (3' 4) 

requires   D   to  lie  inside  the  circle  shown  there. 


A 
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/>- 
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/ 

\    ' 

Figure  3 

The  side  of  the  small  square  equals  the  radius  of  the  circle. 

The  rest  of  this  section  is  about  the  relation  of  the 
stability  region  to  the   C-F-L   condition  for  some  other  dif- 
ference schemes.   Our  starting  point  is  the  result  of  S.  Hahn 
[3]  about  three-point  schemes.   For  these  the   C-P-L   condition 

■X- 

is  sufficient  for  stability  .   Next  we  shall  investigate  four-point 
schemes  of  the  kind  where   S,  v   at   p   is  a  linear  combination 


*  This  result  has  been  extended  in  [9]  to  include  nonsymmetric 
hyperbolic  equations. 
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of  the  values  of   v   at   the  four  lattice  neighbors  of   p: 


(4.4) 


S,  V 

h 


(p)  =  A—   CjV(pj) 


where   p  .   are  the  four  points  Indicated  on  figure  h 

J 


y 


'2 


p. 


p 


'h 


p- 


X 


Figure   4 


We  must  choose  the  coefficients   c  .   so  that  the  difference 
scheme  (4,4)  approximates  the  differential  equation  (l.l) 
at  least  to  first  order.   According  to  Section  2  this 
imposes  three  linear  conditions  on  the  coefficients   c  ; 
it  is  easy  to  show  that  these  linear  relations  are 
satisfied  if  and  only  if   c  .   have  the  form 

J 


(4.5) 


'1,5   ^ 


^  (I  +  M)  t|  ,       c^^^  .  1(1  -  M)  t|  , 


M  an  arbitrary  matrix. 
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For  a  four  point  scheme  as  above  the  limit  of  the 
domains  of  dependence  D   is  the  square  as  shown  in  Figure 
h,    i.e.   with  vertices   (0,±l)   and   (±1,0).   Denote   this 
square  again  by  S;      then  the   C-F-L   condition  is  expressed 
by  (4.2).   As  before,  it  suffices  to  consider  (4.2)  for 
those  values  of   (^,r|)   which  are  perpendicular  to  the 
sides  of   S,  i.e.  (+l,l)   and   (±l,-l).   We  get 


A    (A±B)  <  1 
max^    ^  — 


1   <  A  .  (A+B) 

—  mm   - 


As  before,  these  scalar  inequalities  can  be  expressed  as 
matrix  inequalities: 


(4.6)  (AtB)2  - 


<. 


Introducing  A+B  =  E  and  A-B  =  F  this  can  be  rewritten  as: 


2  .  ^     „2 


(4.6)  E^  <  I,    F^  <  I. 

With  (4.5)  as  choice  for  the  coefficients  the 
amplification  matrix  takes  the  form 

(.4.7;      C  =  — —  cos  fc,  +  —^   cos  r\   +   i(A  sin  |  +  B  sm  ri ) 
=  R  +  i  J  . 
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Assume  now  that  (4.6)  Is  satisfied;  our  aim  is  to 
choose  M   so  that  the  difference  scheme  is  stable.   We 
shall  verify  stability  by  showing  that  for  the  appropriate 
choice  of  M  the  hypothesis  of  theorem  3  is  satisfied. 

To  estimate  (Cu,u)   ,  u  any  unit  vector,  we  write 
as  before 

(4.8)      ( Cu , u )  =  ( Ru , u  )  +  i ( Ju , u  )  =  r  +  i  j , 

By  (4.7)  we  have 
/ 1,  r^\  1+m      *  ,  1-m 

(4,9j        r  =  — ;r—  COS  t,     +    —p—    COS  T] 

where   m  abbreviates 

(4.10)  m  =  (Mu,u)o 

By  the  Schwarz  inequality 

(4.11)  j^  =  (Ju,u)2  <  I  Jul  2  , 

We  rewrite   J,  as  given  by  (4.7)^  in  terms  of   E 
and  P: 


-r    E+P   .   .  ,  E-P   . 

J  =  —;r-    Sm  t,    +  —;r-    Sin  rj , 
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Squaring  we  get 


'■2         „„,„„      ■^„.„2s         _.^2^  \ 


2         ^?  /   sin    e   +    sin   ri  '       ,    EF+FE      '  sin    ^    -    sin   n 
J      =    E     I    p  12  2  / 


,    t:,2  /  sm    g,    -    sm   ri 


Using  inequalities  (4.6)    and  the  abbreviation 


(4.12)         N  =  — p — 


we  get 


o    ^  sin^e  +  sin^ri  ,  .,  sin^^  -  sin  n 
J  <  1  2  2 


Substituting  the  above  into  (4.1l)  we  get 
(4.15)    j'<^  sin2^  +^sin^ 
where   n   abbreviates 


(4,14)    n  =  (Nu,u)„ 

By  (4.8),   |(Cu,u)|^  =  r^  +  j^,   so  to  verify  the 

2    .2 
hypothesis  of  theorem  3  we  have  to  show  that   r  +  j  ^  1^ 

¥e  shall  deal  first  with  the  case  of   ^,'0   small. 


-  38  - 


Expanding  In  Taylor  series  up  to  terms  of  second 
order  we  get  from  (4.9)  and  (4.13) 


-,    1+m  .2        1-m   2 


•2  /  1+n  .2  ,  1-n   2 


plus  terms  of  order  higher  than  two.   Squaring  the  first 
relation  and  adding  It  to  the  second  gives 


2  ,  .2  /  .  ,  n-m  ,  ^2  2\ 


In  order  for  the  right  side  not  to  exceed  one  for  any   ^,rj 
we  must  have   n  =  m   for  all   Uo   In  view  of  the  definitions 
(4.10)  and  (4.14)  this  means  that  we  must  choose   M  =  N. 
By  (4.12)   N  =  M±^  =  A^  -  B^;      with  this  choice   C 
becomes 


(4.7)'   C  -  ^(I  +  A^-b2)cos  I  +  1(1  +  B^-A^)cos  Ti 


+  1(A  sin  i   +  B   sin  rj). 

Theorem  5;   If  the   C-P-L   condition  (4.6)  is  satisfied 
then  the  difference  scheme  associated  with   C   given  by 
(4.7)  is  stable. 
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Proof;  We  shall  prove  that   r^  +  j^  <  1   for  all  real 

i      and   ri.   Squaring  (4.9)  and  adding  it  to  (4.15)  we  get, 
using   n  =  m  : 

(4.15)   r2  -<-  /  ^  (^  '  '  cos^l  +  I±n  sin^e  +   ^  '  cos^ 


1-n  ^,2^  1-n^ 


+  — p—  Sm  'Q   +  -— cos  ^  cos  Tj 


Writing   1  -  cos^l   for   sin^^  etc.,  the  right  side  can  be 
written  as 


1  -  n   /     ^  \  c; 

1  -  — n— —  (cos  ^  -  cos  r\)     . 


So  we  get  the  inequality 


2     2        1-n^  /     ^         n2 
r  +  j   <  1  -  -^—  (cos  i    -    cos  Ti) 


We  claim  that  the  right  side  does  not  exceed  one;  to  see 
this  we  have  to  show  that   |n|   does  not  exceed  one.   Usin^ 
(4,14)  and  (4.12)  we  have 

n  =    (Nu,u)  -  -^^-^   u,u   -  Re(Eu,Fu); 

by  the  Schwarz  inequality  the  absolute  value  of  the  right 
side  does  not  exceed   |Eu[|Fu[   which  does  not  exceed   1 
since   u   is  a  unit  vector  and  by  (4,6)    the  norms  of   E 
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and   F  do  not  exceed  one.   In  view  of  theorem  3  this 

completes  the  proof  of  theorem  5- 

It  would  be  Interesting  to  find  a  more  Intuitive 

t 
way  of  deriving  the  difference  scheme  (4,7)  . 

It  Is  not  clear  whether  the  greater  stability  and 

simplicity  of  the  four  point  scheme  (4.7)  compensates  for 

the  loss  of  accuracy  as  compared  to  the  schemes  discussed 

in  Section  3. 
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5.    Stability  of  Difference  Schemes^  with  Variable  Coef- 
ficients 

Since  the  use  of  difference  schemes  is  for  solving 
differential  equations  with  variable  coefficients  (and 
nonlinear  equations)  it  is  essential  to  prove  that  the 
difference  schemes  devised  in  section  3  are  stable  also 
when  the  coefficients  are  variable.   Experience  both 
practical  and  theoretical  shows  that  the  stability  of  dif- 
ference schemes  with  smoothly  varying  coefficients  is 

governed  by  the  local  amipllf ication  matrix   C  (x,|) 

o 

defined  as 


C  (x,0  =  XI  c  (x,0)e^J^ 
o  .         J 


where   c.(x,0)   denotes  the  value  of   c.   for   h  =  0. 
J  J 

A  theorem  in  this  direction  is  the  following,  see  [  12]: 

Theorem  6 ;   If  for  every  x  and   C   the  norm  of  the 
local  amplification  matrix  does  not  exceed  one  then  the 
associated  difference  scheme  is  stable. 

Actually,  as  shown  in  [lOj,a  slightly  stronger  result 
holds;  it  is  sufficient  to  require  that  the  hypothesis 
hold  for  i,r\      sufficiently  small  if  in  addition  the 
eigenvalues  of   C   are  less  than  one  for   |^|  +  \r\\ 
bounded   away  from  zero. 

It  is  not  hard  to  show  that  the  first  difference 
schemie  discussed  in  Section  3  does  not  satisfy  the 


4^ 


hypothesis  of  theorem  6,    even  when  weakened  as  above, 
unless  the  coeiTlclent  matrices  A   and  3   satisfy  some 
rather  restrictive  conditions  such  as  being  bounded  from 
below  as  well  as  from  above,  or  commuting  with  each  other, 
The  second  difference  scheme  on  the  other  hand  does  not 
satisfy  the  hypothesis  of  theorem  6,  at  least  If  the 
stability  condition  {j>A)    Is  replaced  by  a  slightly 
stronger  one. 

We  surmise  that  the  following  result  holds: 

Conjecture ;   if  for  every  x  and   ^   tne  local 
amplification  matrix  satisfies  the  hypothesis  of  theorem 
3  in  Section  2,  then  the  associated  difference  scheme  Is 
stable . 
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6.    Systems  of  Conservation  laws. 

In  this  section  we  shall  adapt  the  difference  schemes 
described  in  Section  3  to  the  construction  of  approxlm.ate 
weak  solutions*  of  systems  of  conservation  laws,  i.e.  of 
equations  of  the  form 


'6.1)  u,  =  f  +  g  , 


u   a  vector  of  unknown  functions,   f   and   g  nonlinear 

vector  valued  functions  of   u.   Carrying  out  the  differentiation 

on  the  right  brings  (6.1)  into  the  form 

(6.1)  u,  =  Au   +  Bu 

t     X     y 

where 

(6.2)  A  -  grad  f,      B  =  grad  g  . 

We  assume  that  the  matrices   A   and   B  can  be  made  symmetric 
by  the  same  similarity  transformation;  this  guarantees  that 
(6.1)   is  hyperbolic. 

It  was  shown  in  [ll],  for  systems  of  conservation  laws 
in  one  space  variable  -  and  the  proof  carries  over  to  several 
variables  -  that  if  we  approximate  (6.l)  by  a  difference 
equation  in  conservation  form,  then  the  strong  limit  of  the 
approximate  solution  is  a  weak  solution  of  the  conservation 


* 

See  e.g. [8]  for  a  discussion  of  the  theory  of  weak  solutions 
of  systems  of  conservation  laws. 


-  kh  - 


law.   By  conservation  form  we  mean 


(6.3)      v(t+h)  -  v(t)  +  D^F  +  D^G   , 


h        h 
where   D   and   D   denote  the  centered  difference  operators 


(6.4)      (d\)(x)  =  u(x  +  h/2)  -  u(x  -  h/2) 


and  similarly  In   y,   and  where   P  and   G  denote  functions  of 
the  values  of  T^  ^^u,   j   ranging  over  some  finite  set  such  that 
If  all  the  arguments   T^  ^u  are  set  equal;   P  reduces  to  f, 
G  to   g.   In  [ll]  we  have  shown  how  to  construct  In  the  case 
of  one  space  variable  difference  equations  In  conservation  form 
which  approximate  a  given  system  of  conservation  law  with 

second  order  accuracy.   Here  we  extend  this  to  two  space 
variables . 

Following  the  method  described  In  Section  1  we  write: 

(6.5)     u(t+h)  ^  u  +  hu,  +  -^  u^, 

modulo  terms  of  third  order  In  h.   Differentiating  (60 l)  with 
respect  to  t  we  get^  using  (6.2) 


'\t   ^^tx  +  Sty  =  ^^'\^x   +  (^-^t) 


y 


-  (Af   +  Ag  )   +  (Ef   +  Bg  ) 

^  X    y  X   ^  X    y  y 
Substituting  this  and  (6.I)  into  (6.5)  gives 
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,2 

(6.6)      u(t+h)  -  u  +  [hf  +  ^  (Af   +  Ag  )] 

d  ^     y   X 


.2 

+  [hg  +  ^  (Bf   +  Bg  )] 


The  important  point  about  the  above  formula  is  that  the  right 
side  is  in  conservation  form  and  so  can  be  approximated  by  a 
difference  expression  of  the  same  kind.   Using  the  abbreviation 


(6.7)      m\i  =  [u(x  +  h/2)  +  u(x  -  h/2)/2,  etc 


we   approximate   the   right    side   of    (6,6)   by      S  u     defined   as 
follows : 


(6.8)  S^^u  =   u  +   dVt   +   oVg   +   \  D^^AD^f 

h  XX  y   y  "^  ■   X      X 


+    \  D^^-BD^g   +    i   D'^M^^AD'^M^g   +   \   D^M^BD^M^V, 
'^yy  2xxyy'='         2yyxx 


Clearly  with      S        defined   by    (6.8), 


;6.9)  v(t   +   h)    =   S^v 


is  of  the  general  form  (6  =,5). 

It  is  easy  to  see  that  (6.9)  is  accurate  to  second  order, 
and  that  the  amplification  matrix  associated  with  the  linearized 
form  of  (6,9)  is  given  by  (3,3)-   This  indicates  that  (6.9)  is 
stable,  at  least  away  from  regions  where   u   is  discontinuous. 

-  46  - 


T 

Define   S,   as 
h 


(6.8)'     S '  =  S^  -  i  D'^(D^-)^AD^r  -  i  D^iD^fBB^-g    ; 


the  amplification  matrix  associated  with  the  linearized  form 
of 


(6.9)'    v(t+h)  =  S,V 


is  (3o)  •   (6.9)  and  (6.9)    are  our  proposed  difference 
schemes . 

The  main  use  of  difference  equations  of  the  form  (6.9) 
and  (6.9)    is  to  calculate  approximations  to  discontinuous 
solutions  of  (6.1),   In  such  calculations  it  is  Important  to 
keep  discontinuities  In  the  approximate  solutions  fairly  sharp. 
The  method  of  artificial  viscosity  [13 ]  was  developed  to 
accomplish  this,  and  in  [ll]  we  have  devised  a  way  of 
introducing  artificial  viscosity  for  arbitrary  systems  of 
conservation  laws  in  one  space  variable.   We  hope  in  a  future 
publication  to  do  the  same  in  the  case  of  two  space  variables. 
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